# author: Jochen Rehmert, University of Basel
# journal: Journal of Politics
# article: Intra-party competition, geographic responsiveness 
#           and incumbent deselection in closed-list PR
# date: 12th August 2024
# content: Script to replicate Figures 4, 5, 13 and 
# Tables 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 
# in the main text and the Appendix
# needs data: "obs.RDS"
# see also: "data_description_obs.txt"

# directory
setwd("")

sink("log_rep_tables4_6-15_figures4_5_13.txt")
cat("start logging\n")

# data & packages
library(lmtest)
library(multiwayvcov)
library(sandwich)
library(stargazer)
library(clubSandwich)

dat <- readRDS("obs.rds")


#create function for calculation of observed-values approach

ame.p <- function(input, dv = c("rank_diff", "ranks_competed", "margin_final"),
                  iv = c("pq_both_coverage_land","pq_both_num_land"), 
                  seed = 1104, nsim = 1000){
  require(sandwich);require(lmtest);require(MASS)
  dat.tmp <- input[["model"]]
  colnames(dat.tmp)[which(colnames(dat.tmp) %in% "as.factor(incumbent)")] <- "incumbent"
  if(dv %in% c("rank_diff","margin_final") ){form1 <- as.formula(paste(dv,"~", paste(names(dat.tmp[-1]), collapse ="+") ))}else{form1 <- input[["formula"]]}
  #if(dv %in% c("rank_final_prob") ){form1 <- as.formula(paste(dv,"~", paste(names(dat.tmp[-1])[-which(names(dat.tmp[-1]) == "smd_nomination")], collapse ="+") ))}else{form1 <- input[["formula"]]}
  
  # run models
  if(dv == "f_denied"){mod1 <- glm(form1, data = dat.tmp, family = binomial("logit"))}
  #if(dv == "rank_final_prob"){mod1 <- glm(form1, data = dat.tmp, family = quasibinomial("logit"))}
  if(dv %in% c("rank_diff","margin_final")){mod1 <- lm(form1, data = dat.tmp)}
  if(dv == "ranks_competed"|dv == "numcand_final"){mod1 <- glm(form1, data = dat.tmp, family = poisson("log"))}
  
  require(fastDummies)
  mod.dat.0 <- mod.dat.1 <- data.frame(cbind(dummy_cols(dat.tmp$incumbent)[-1], mod1[["model"]][-c(1:2)]))
  mod.dat.1[[iv]] <- dat.tmp[[iv]] + mean(dat.tmp[[iv]], na.rm = TRUE)
  
  # simulations mod1 cov(mod1)
  set.seed(seed)
  
  b.sim1 <- mvrnorm(nsim, coef(mod1), vcovHC(mod1, type="HC"))
 
  lin.preds.0 <- as.matrix(mod.dat.0) %*% t(b.sim1)  
  lin.preds.1 <- as.matrix(mod.dat.1) %*% t(b.sim1)  
  
  if(dv %in% c("f_denied")){   # , "rank_final_prob"
    preds.0 <- exp(lin.preds.0)/(1+exp(lin.preds.0))
    preds.1 <- exp(lin.preds.1)/(1+exp(lin.preds.1))
  }
  if(dv %in% c("ranks_competed","numcand_final")){
    preds.0 <- exp(lin.preds.0)
    preds.1 <- exp(lin.preds.1)
  }
  if(dv %in% c("rank_diff","margin_final")){
    preds.0 <- lin.preds.0
    preds.1 <- lin.preds.1
  }
  
  preds.0.mean <- apply(preds.0, 2, mean, na.rm = TRUE)
  preds.1.mean <- apply(preds.1, 2, mean, na.rm = TRUE)
  
  ame.95 <- quantile(preds.1.mean - preds.0.mean, probs = c(0.5, 0.025, 0.975), na.rm = TRUE)
  
  ame.90 <- quantile(preds.1.mean - preds.0.mean, probs = c(0.5, 0.05, 0.95), na.rm = TRUE)
  
  return(list(iv = iv,
              dv = dv,
              ame90 = ame.90,
              ame95 = ame.95))
  
}



#################################################
# 1) number of ranks competed
#################################################

library(AER)
##### -- All -- #####
### --- PQs - Share --- ###
summary(a2.c <- glm(ranks_competed ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                      office + elected+ theta_mean_land_deviation + ln_pq_total+ pq_both_coverage_land ,  data = dat, family = poisson("log")))
a2.c.se.rob <- coeftest(a2.c, vcov = vcovHC(a2.c, "HC"))[,2];coeftest(a2.c, vcov = vcovHC(a2.c, "HC"))
a2.c.se.cl <- coef_test(a2.c, vcov = "CR2", 
                        cluster = a2.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
rank.cov.a <- ame.p(a2.c, dv = "ranks_competed", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(a2.d <- glm(ranks_competed ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                      office + elected+ theta_mean_land_deviation + ln_pq_total+ ln_pq_land ,  data = dat, family = poisson("log")))
a2.d.se.rob <- coeftest(a2.d, vcov = vcovHC(a2.d, "HC"))[,2];coeftest(a2.d, vcov = vcovHC(a2.d, "HC"))
a2.d.se.cl <- coef_test(a2.d, vcov = "CR2", 
                        cluster = a2.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
rank.num.a <- ame.p(a2.d, dv = "ranks_competed", iv = "ln_pq_land",seed = 1104, nsim = 1000)


##### -- Women -- #####
### --- PQs - Share --- ###
# remove office
summary(w2.c <- glm(ranks_competed ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                       elected+ theta_mean_land_deviation + ln_pq_total+ pq_both_coverage_land ,  data = dat[dat$sex == "female",], family = poisson("log")))
w2.c.se.rob <- coeftest(w2.c, vcov = vcovHC(w2.c, "HC"))[,2];coeftest(w2.c, vcov = vcovHC(w2.c, "HC"))
w2.c.se.cl <- coef_test(w2.c, vcov = "CR2", 
          cluster = w2.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
rank.cov.w <- ame.p(w2.c, dv = "ranks_competed", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(w2.d <- glm(ranks_competed ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                       elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land,  data = dat[dat$sex == "female",], family = poisson("log")))
w2.d.se.rob <- coeftest(w2.d, vcov = vcovHC(w2.d, "HC"))[,2];coeftest(w2.d, vcov = vcovHC(w2.d, "HC"))
w2.d.se.cl <- coef_test(w2.d, vcov = "CR2", 
                            cluster = w2.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
rank.num.w <- ame.p(w2.d, dv = "ranks_competed", iv = "ln_pq_land",seed = 1104, nsim = 1000)



##### -- Men -- #####
### --- PQs - Share --- ###
# remove smd_nomination
summary(m2.c <- glm(ranks_competed ~ - 1+ as.factor(incumbent) + rank_before_prob +  
                      office + elected+ theta_mean_land_deviation+ ln_pq_total +pq_both_coverage_land ,  data = dat[dat$sex == "male",], family = poisson("log")))
m2.c.se.rob <- coeftest(m2.c, vcov = vcovHC(m2.c, "HC"))[,2];coeftest(m2.c, vcov = vcovHC(m2.c, "HC"))
m2.c.se.cl <- coef_test(m2.c, vcov = "CR2", 
                        cluster = m2.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
rank.cov.m <- ame.p(m2.c, dv = "ranks_competed", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(m2.d <- glm(ranks_competed ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                      office + elected+ theta_mean_land_deviation +ln_pq_total + ln_pq_land,  data = dat[dat$sex == "male",], family = poisson("log")))
m2.d.se.rob <- coeftest(m2.d, vcov = vcovHC(m2.d, "HC"))[,2];coeftest(m2.d, vcov = vcovHC(m2.d, "HC"))
m2.d.se.cl <- coef_test(m2.d, vcov = "CR2", 
                        cluster = m2.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
rank.num.m <- ame.p(m2.d, dv = "ranks_competed", iv = "ln_pq_land",seed = 1104, nsim = 1000)


################
# Table 6 in Appendix
################
# all incumbents
cat("create table 6\n")
stargazer(list( a2.c, a2.d, a2.c, a2.d),
          se = list( a2.c.se.rob, a2.d.se.rob, 
                     a2.c.se.cl, a2.d.se.cl), omit = c("incumbent"),
          covariate.labels = c("Pr (Election) Rank t-1",
                               "Office Holder",
                               "Times Elected",
                               "Ideological Deviation from Landesverband",
                               "Total Number of PQs, log",
                               "Share of Kreise Covered",
                               "Number of Land Locations Mentioned in PQs, log"))

##########££####
# Table 9 in Appendix
################
# split by gender
cat("create table 9\n")
stargazer(list( w2.c, m2.c, w2.d, m2.d,w2.c, m2.c, w2.d, m2.d),
          se = list( w2.c.se.rob, m2.c.se.rob, w2.d.se.rob, m2.d.se.rob,
                     w2.c.se.cl, m2.c.se.cl, w2.d.se.cl, m2.d.se.cl), 
          omit = c("incumbent"),
          covariate.labels = c("Pr (Election) Rank t-1",
                               "Office Holder",
                               "Times Elected",
                               "Ideological Deviation from Landesverband",
                               "Total Number of PQs, log",
                               "Share of Kreise Covered",
                               "Number of Land Locations Mentioned in PQs, log"))

# test for overdispersion
dispersiontest(a2.c)
dispersiontest(a2.d)
dispersiontest(w2.c)
dispersiontest(w2.d)
dispersiontest(m2.c)
dispersiontest(m2.d)
#################################################
# 2) difference in rank
#################################################
dat$rank_diff <- dat$rank_before - dat$rank_final

dat$rank_diff[!is.na(dat$rank_diff) &  dat$rank_final == 0] <- -1

##### -- All -- #####
### --- PQs - Share --- ###
summary(a4.c <- lm(rank_diff ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation+ ln_pq_total + pq_both_coverage_land ,  data = dat))
a4.c.se.rob <- coeftest(a4.c, vcov = vcovHC(a4.c, "HC"))[,2];coeftest(a4.c, vcov = vcovHC(a4.c, "HC"))
a4.c.se.cl <- coef_test(a4.c, vcov = "CR2", 
                        cluster = a4.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
diff.cov.a <- ame.p(a4.c, dv = "rank_diff", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(a4.d <- lm(rank_diff ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land  ,  data = dat))
a4.d.se.rob <- coeftest(a4.d, vcov = vcovHC(a4.d, "HC"))[,2];coeftest(a4.d, vcov = vcovHC(a4.d, "HC"))
a4.d.se.cl <- coef_test(a4.d, vcov = "CR2", 
                        cluster = a4.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
diff.num.a <- ame.p(a4.d, dv = "rank_diff", iv = "ln_pq_land",seed = 1104, nsim = 1000)



##### -- Women -- #####
### --- PQs - Share --- ###
summary(w4.c <- lm(rank_diff ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation+ ln_pq_total + pq_both_coverage_land ,  data = dat[dat$sex == "female",]))
w4.c.se.rob <- coeftest(w4.c, vcov = vcovHC(w4.c, "HC"))[,2];coeftest(w4.c, vcov = vcovHC(w4.c, "HC"))
w4.c.se.cl <- coef_test(w4.c, vcov = "CR2", 
          cluster = w4.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
diff.cov.w <- ame.p(w4.c, dv = "rank_diff", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(w4.d <- lm(rank_diff ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land  ,  data = dat[dat$sex == "female",]))
w4.d.se.rob <- coeftest(w4.d, vcov = vcovHC(w4.d, "HC"))[,2];coeftest(w4.d, vcov = vcovHC(w4.d, "HC"))
w4.d.se.cl <- coef_test(w4.d, vcov = "CR2", 
                        cluster = w4.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
diff.num.w <- ame.p(w4.d, dv = "rank_diff", iv = "ln_pq_land",seed = 1104, nsim = 1000)


##### -- Men -- #####
### --- PQs - Share --- ###
summary(m4.c <- lm(rank_diff ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total + pq_both_coverage_land ,  data = dat[dat$sex == "male",]))#, family = binomial("logit")))
m4.c.se.rob <- coeftest(m4.c, vcov = vcovHC(m4.c, "HC"))[,2];coeftest(m4.c, vcov = vcovHC(m4.c, "HC"))
m4.c.se.cl <- coef_test(m4.c, vcov = "CR2", 
          cluster = m4.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
diff.cov.m <- ame.p(m4.c, dv = "rank_diff", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(m4.d <- lm(rank_diff ~ - 1+ as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation+ ln_pq_total + ln_pq_land  ,  data = dat[dat$sex == "male",]))#, family = binomial("logit")))
m4.d.se.rob <- coeftest(m4.d, vcov = vcovHC(m4.d, "HC"))[,2];coeftest(m4.d, vcov = vcovHC(m4.d, "HC"))
m4.d.se.cl <- coef_test(m4.d, vcov = "CR2", 
          cluster = m4.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
diff.num.m <- ame.p(m4.d, dv = "rank_diff", iv = "ln_pq_land",seed = 1104, nsim = 1000)


################
# Table 7 in Appendix
################
# all
cat("create table 7\n")
stargazer(list( a4.c, a4.d,  a4.c, a4.d),
          se = list( a4.c.se.rob, a4.d.se.rob, 
                     a4.c.se.cl, a4.d.se.cl), omit = c("incumbent")
          ,
          covariate.labels = c("Pr (Election) Rank t-1",
                               "Office Holder",
                               "Times Elected",
                               "Ideological Deviation from Landesverband",
                               "Total Number of PQs, log",
                               "Share of Kreise Covered",
                               "Number of Land Locations Mentioned in PQs, log"))

################
# Table 10 in Appendix
################
# split by gender
cat("create table 10\n")
stargazer(list( w4.c, m4.c, w4.d, m4.d,w4.c, m4.c, w4.d, m4.d),
          se = list( w4.c.se.rob, m4.c.se.rob, w4.d.se.rob, m4.d.se.rob,
                     w4.c.se.cl, m4.c.se.cl, w4.d.se.cl, m4.d.se.cl), 
          omit = c("incumbent"),
          covariate.labels = c("Pr (Election) Rank t-1",
                               "Office Holder",
                               "Times Elected",
                               "Ideological Deviation from Landesverband",
                               "Total Number of PQs, log",
                               "Share of Kreise Covered",
                               "Number of Land Locations Mentioned in PQs, log"))


#################################################
# 3) margin for final rank
#################################################

##### -- All -- #####
### --- PQs - Share --- ###
summary(a1.c <- lm(margin_final ~ -1 +  as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total + pq_both_coverage_land ,  data =  dat))
a1.c.se.rob <- coeftest(a1.c, vcov = vcovHC(a1.c, "HC"))[,2];coeftest(a1.c, vcov = vcovHC(a1.c, "HC"))
a1.c.se.cl <- coef_test(a1.c, vcov = "CR2", 
                        cluster = a1.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
mar.cov.a <- ame.p(a1.c, dv = "margin_final", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(a1.d <- lm(margin_final ~ -1 +  as.factor(incumbent) + rank_before_prob +
                     office + elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land    ,  data = dat))
a1.d.se.rob <- coeftest(a1.d, vcov = vcovHC(a1.d, "HC"))[,2];coeftest(a1.d, vcov = vcovHC(a1.d, "HC"))
a1.d.se.cl <- coef_test(a1.d, vcov = "CR2", 
                        cluster = a1.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
mar.num.a <- ame.p(a1.d, dv = "margin_final", iv = "ln_pq_land",seed = 1104, nsim = 1000)


##### -- Women -- #####
### --- PQs - Share --- ###
summary(w1.c <- lm(margin_final ~ -1 +  as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total + pq_both_coverage_land  ,  data =  dat[dat$sex == "female",]))
w1.c.se.rob <- coeftest(w1.c, vcov = vcovHC(w1.c, "HC"))[,2];coeftest(w1.c, vcov = vcovHC(w1.c, "HC"))
w1.c.se.cl <- coef_test(w1.c, vcov = "CR2", 
          cluster = w1.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
mar.cov.w <- ame.p(w1.c, dv = "margin_final", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(w1.d <- lm(margin_final ~ -1 +  as.factor(incumbent) + rank_before_prob +
                     office + elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land   ,  data = dat[dat$sex == "female",]))
w1.d.se.rob <- coeftest(w1.d, vcov = vcovHC(w1.d, "HC"))[,2];coeftest(w1.d, vcov = vcovHC(w1.d, "HC"))
w1.d.se.cl <- coef_test(w1.d, vcov = "CR2", 
          cluster = w1.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
mar.num.w <- ame.p(w1.d, dv = "margin_final", iv = "ln_pq_land",seed = 1104, nsim = 1000)


##### -- Men -- #####
### --- PQs - Share --- ###
summary(m1.c <- lm(margin_final ~ -1 +  as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total+ pq_both_coverage_land   ,  data =  dat[dat$sex == "male",]))
m1.c.se.rob <- coeftest(m1.c, vcov = vcovHC(m1.c, "HC"))[,2];coeftest(m1.c, vcov = vcovHC(m1.c, "HC"))
m1.c.se.cl <- coef_test(m1.c, vcov = "CR2", 
          cluster = m1.c[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
mar.cov.m <- ame.p(m1.c, dv = "margin_final", iv = "pq_both_coverage_land",seed = 1104, nsim = 1000)

### --- PQs - Number --- ###
summary(m1.d <- lm(margin_final ~ -1 +  as.factor(incumbent) + rank_before_prob +
                     office + elected+ ln_pq_total + ln_pq_land + theta_mean_land_deviation,  data = dat[dat$sex == "male",]))
m1.d.se.cl <- coef_test(m1.d, vcov = "CR2", 
          cluster = m1.d[["model"]]$`as.factor(incumbent)`, test = "Satterthwaite")[,3]
m1.d.se.rob <- coeftest(m1.d, vcov = vcovHC(m1.d, "HC"))[,2];coeftest(m1.d, vcov = vcovHC(m1.d, "HC"))

mar.num.m <- ame.p(m1.d, dv = "margin_final", iv = "ln_pq_land",seed = 1104, nsim = 1000)

################
# Table 8 in Appendix
################
# all
cat("create table 8\n")
stargazer(list(  a1.c, a1.d,  a1.c, a1.d),
          se = list(  a1.c.se.rob, a1.d.se.rob, 
                      a1.c.se.cl, a1.d.se.cl), 
          omit = c("incumbent"),
          covariate.labels = c("Pr (Election) Rank t-1",
                               "Office Holder",
                               "Times Elected",
                               "Ideological Deviation from Landesverband",
                               "Total Number of PQs, log",
                               "Share of Kreise Covered",
                               "Number of Land Locations Mentioned in PQs, log"))

################
# Table 11 in Appendix
################
# split by gender
cat("create table 11\n")
stargazer(list(  w1.c, m1.c, w1.d, m1.d,w1.c, m1.c, w1.d, m1.d),
          se = list(  w1.c.se.rob, m1.c.se.rob, w1.d.se.rob, m1.d.se.rob,
                      w1.c.se.cl, m1.c.se.cl, w1.d.se.cl, m1.d.se.cl),
          omit = c("incumbent"),
          covariate.labels = c("Pr (Election) Rank t-1",
                               "Office Holder",
                               "Times Elected",
                               "Ideological Deviation from Landesverband",
                               "Total Number of PQs, log",
                               "Share of Kreise Covered",
                               "Number of Land Locations Mentioned in PQs, log"))

#################################################
#  election probability of final rank
#################################################

##### -- All -- #####
### --- PQs - Share --- ###
summary(a1.c <- glm(rank_final_prob ~ -1 +as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total + pq_both_coverage_land  , family = quasibinomial("logit"),  data =  dat))
a1.c.se.rob <- coeftest(a1.c, vcov = vcovHC(a1.c, "HC"))[,2];coeftest(a1.c, vcov = vcovHC(a1.c, "HC"))
a1.c.se.cl <- coef_test(a1.c, vcov = "CR2",  cluster = dat$land, test = "Satterthwaite")[,3]


### --- PQs - Number --- ###
summary(a1.d <- glm(rank_final_prob ~ -1 +as.factor(incumbent) + rank_before_prob +
                     office + elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land   , family = quasibinomial("logit"), data = dat))
a1.d.se.rob <- coeftest(a1.d, vcov = vcovHC(a1.d, "HC"))[,2];coeftest(a1.d, vcov = vcovHC(a1.d, "HC"))
a1.d.se.cl <- coef_test(a1.d, vcov = "CR2",   cluster = dat$land, test = "Satterthwaite")[,3]


##### -- Women -- #####
### --- PQs - Share --- ###
summary(w1.c <- glm(rank_final_prob ~ -1 +as.factor(incumbent) + rank_before_prob +
                     office + elected+ theta_mean_land_deviation + ln_pq_total + pq_both_coverage_land  , family = quasibinomial("logit"), data =  dat[dat$sex == "female",]))
w1.c.se.rob <- coeftest(w1.c, vcov = vcovHC(w1.c, "HC"))[,2];coeftest(w1.c, vcov = vcovHC(w1.c, "HC"))
w1.c.se.cl <- coef_test(w1.c, vcov = "CR2", cluster = dat$land[dat$sex == "female"], test = "Satterthwaite")[,3]
 

### --- PQs - Number --- ###
summary(w1.d <- glm(rank_final_prob ~ -1 +as.factor(incumbent)+ rank_before_prob +
                     office + elected+ theta_mean_land_deviation + ln_pq_total + ln_pq_land   , family = quasibinomial("logit"),  data = dat[dat$sex == "female",]))
w1.d.se.rob <- coeftest(w1.d, vcov = vcovHC(w1.d, "HC"))[,2];coeftest(w1.d, vcov = vcovHC(w1.d, "HC"))
w1.d.se.cl <- coef_test(w1.d, vcov = "CR2",  cluster = dat$land[dat$sex == "female"], test = "Satterthwaite")[,3]

 
##### -- Men -- #####
### --- PQs - Share --- ###
summary(m1.c <- glm(rank_final_prob ~ -1 +as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total+ pq_both_coverage_land   , family = quasibinomial("logit"), data =  dat[dat$sex == "male",]))
m1.c.se.rob <- coeftest(m1.c, vcov = vcovHC(m1.c, "HC"))[,2];coeftest(m1.c, vcov = vcovHC(m1.c, "HC"))
m1.c.se.cl <- coef_test(m1.c, vcov = "CR2",  cluster = dat$land[dat$sex == "male"], test = "Satterthwaite")[,3]

 
### --- PQs - Number --- ###
summary(m1.d <- glm(rank_final_prob ~ -1 +as.factor(incumbent) + rank_before_prob +
                     office + elected+ ln_pq_total + ln_pq_land + theta_mean_land_deviation, family = quasibinomial("logit"),  data = dat[dat$sex == "male",]))
m1.d.se.cl <- coef_test(m1.d, vcov = "CR2",  cluster = dat$land[dat$sex == "male"], test = "Satterthwaite")[,3]
m1.d.se.rob <- coeftest(m1.d, vcov = vcovHC(m1.d, "HC"))[,2];coeftest(m1.d, vcov = vcovHC(m1.d, "HC"))
 

################
# Table 13 in Appendix
################
# split by gender
cat("create table 13\n")
stargazer(list(  w1.c, m1.c, w1.d, m1.d,w1.c, m1.c, w1.d, m1.d),
          se = list(  w1.c.se.rob, m1.c.se.rob, w1.d.se.rob, m1.d.se.rob,
                      w1.c.se.cl, m1.c.se.cl, w1.d.se.cl, m1.d.se.cl), omit = "incumbent",
          covariate.labels = c("Pr (Election) Rank t-1",
                               "Office Holder",
                               "Times Elected",
                               "Ideological Deviation from Landesverband",
                               "Total Number of PQs, log",
                               "Share of Kreise Covered",
                               "Number of Land Locations Mentioned in PQs, log"),
         add.lines = list(c("Deviance", round(w1.c$deviance,2),
                             round(m1.c$deviance,2), round(w1.d$deviance,2), 
                             round(m1.d$deviance,2),round(w1.c$deviance,2), 
                             round(m1.c$deviance,2), round(w1.d$deviance,2), round(m1.d$deviance,2) )))

#################################################
#                 Descriptives                  #
#################################################
dat.desc <- dat[which(!is.na(dat$theta_mean_land_deviation) & !is.na(dat$rank_diff)),
                c("rank_diff","margin_final", "ranks_competed",  
                  "pq_both_coverage_land","pq_both_num_land" , "pq_total",  
                  "office","elected","theta_mean_land_deviation","rank_before_prob")]

################
# Table 4 in Appendix
################
cat("create table 4\n")
stargazer(dat.desc)



#################################################
#                 plotting AMEs                 #
#################################################

pData <- data.frame(rbind(cbind(outcome = "No. of Ranks Competed"  , iv = "Share of \nKreise Covered", ame = rank.cov.a$ame95[1], ci = 95, lower = rank.cov.a$ame95[2], upper = rank.cov.a$ame95[3] ),
                          cbind(outcome = "No. of Ranks Competed"  , iv = "Share of \nKreise Covered", ame = rank.cov.a$ame90[1], ci = 90, lower = rank.cov.a$ame90[2], upper = rank.cov.a$ame90[3] ),
                          cbind(outcome = "No. of Ranks Competed"  , iv = "Number of \nLand Locations", ame = rank.num.a$ame95[1], ci = 95, lower = rank.num.a$ame95[2], upper = rank.num.a$ame95[3] ),
                          cbind(outcome = "No. of Ranks Competed"  , iv = "Number of \nLand Locations", ame = rank.num.a$ame90[1], ci = 90, lower = rank.num.a$ame90[2], upper = rank.num.a$ame90[3] ),
                          cbind(outcome = "Difference in Rank Position"  , iv = "Share of \nKreise Covered", ame = diff.cov.a$ame95[1], ci = 95, lower = diff.cov.a$ame95[2], upper = diff.cov.a$ame95[3] ),
                          cbind(outcome = "Difference in Rank Position"  , iv = "Share of \nKreise Covered", ame = diff.cov.a$ame90[1], ci = 90, lower = diff.cov.a$ame90[2], upper = diff.cov.a$ame90[3] ),
                          cbind(outcome = "Difference in Rank Position"  , iv = "Number of \nLand Locations", ame = diff.num.a$ame95[1], ci = 95, lower = diff.num.a$ame95[2], upper = diff.num.a$ame95[3] ),
                          cbind(outcome = "Difference in Rank Position"  , iv = "Number of \nLand Locations", ame = diff.num.a$ame90[1], ci = 90, lower = diff.num.a$ame90[2], upper = diff.num.a$ame90[3] ),
                          cbind(outcome = "Vote Margin for Final Rank"  , iv = "Share of \nKreise Covered", ame = mar.cov.a$ame95[1], ci = 95, lower = mar.cov.a$ame95[2], upper = mar.cov.a$ame95[3] ),
                          cbind(outcome = "Vote Margin for Final Rank"  , iv = "Share of \nKreise Covered", ame = mar.cov.a$ame90[1], ci = 90, lower = mar.cov.a$ame90[2], upper = mar.cov.a$ame90[3] ),
                          cbind(outcome = "Vote Margin for Final Rank"  , iv = "Number of \nLand Locations", ame = mar.num.a$ame95[1], ci = 95, lower = mar.num.a$ame95[2], upper = mar.num.a$ame95[3] ),
                          cbind(outcome = "Vote Margin for Final Rank"  , iv = "Number of \nLand Locations", ame = mar.num.a$ame90[1], ci = 90, lower = mar.num.a$ame90[2], upper = mar.num.a$ame90[3] )))
                          

row.names(pData) <- 1:nrow(pData)
pData$ame <- as.numeric(as.character(pData$ame))
pData$lower <- as.numeric(as.character(pData$lower))
pData$upper <- as.numeric(as.character(pData$upper))
pData$ci <- as.numeric(as.character(pData$ci))
pData$outcome <- as.character(pData$outcome)
pData$outcome_f <- factor(pData$outcome, levels = c("No. of Ranks Competed","Vote Margin for Final Rank","Difference in Rank Position"))


library(ggplot2)

p = ggplot(pData[pData$ci == 95 ,])
p = p + geom_pointrange(data = pData[pData$ci == 90 ,], aes(x = iv, y = ame, ymin = lower, ymax = upper))
p = p + geom_pointrange(data = pData[pData$ci == 95 ,],aes(x = iv, y = ame, ymin = lower, ymax = upper), lty = 2)
p = p + geom_hline(aes(yintercept = 0), linetype = 2, lwd = 1.05) + xlab("") + ylab("AME")
p = p + facet_wrap(vars(outcome_f), scales = 'free_x')
p = p + coord_flip() + theme(axis.text = element_text(size = rel(1), face = "bold", colour = "black"),
                             axis.ticks = element_line(colour = "black"),
                             legend.key = element_rect(fill = "white"),
                             plot.margin = unit(c(0, 0, 0, 0), "cm"),
                             panel.background = element_rect(fill = "white", colour = NA),
                             panel.border = element_rect(fill = NA, colour = "grey50"),
                             strip.background = element_rect(fill = "white", colour = "grey50"),
                             strip.text.x = element_text(size = 12, colour = "black", face = "bold"),
                             strip.text.y = element_text(size = 10, colour = "black", face = "bold"))  
################
# Figure 4 in main text
################
cat("create figure 4\n")
p




pData <- data.frame(rbind(cbind(gender = "Woman",outcome = "No. of Ranks Competed"  , iv = "Share of \nKreise Covered", ame = rank.cov.w$ame95[1], ci = 95, lower = rank.cov.w$ame95[2], upper = rank.cov.w$ame95[3] ),
                          cbind(gender = "Woman",outcome = "No. of Ranks Competed"  , iv = "Share of \nKreise Covered", ame = rank.cov.w$ame90[1], ci = 90, lower = rank.cov.w$ame90[2], upper = rank.cov.w$ame90[3] ),
                          cbind(gender = "Woman",outcome = "No. of Ranks Competed"  , iv = "Number of \nLand Locations", ame = rank.num.w$ame95[1], ci = 95, lower = rank.num.w$ame95[2], upper = rank.num.w$ame95[3] ),
                          cbind(gender = "Woman",outcome = "No. of Ranks Competed"  , iv = "Number of \nLand Locations", ame = rank.num.w$ame90[1], ci = 90, lower = rank.num.w$ame90[2], upper = rank.num.w$ame90[3] ),
                          cbind(gender = "Woman",outcome = "Difference in Rank Position"  , iv = "Share of \nKreise Covered", ame = diff.cov.w$ame95[1], ci = 95, lower = diff.cov.w$ame95[2], upper = diff.cov.w$ame95[3] ),
                          cbind(gender = "Woman",outcome = "Difference in Rank Position"  , iv = "Share of \nKreise Covered", ame = diff.cov.w$ame90[1], ci = 90, lower = diff.cov.w$ame90[2], upper = diff.cov.w$ame90[3] ),
                          cbind(gender = "Woman",outcome = "Difference in Rank Position"  , iv = "Number of \nLand Locations", ame = diff.num.w$ame95[1], ci = 95, lower = diff.num.w$ame95[2], upper = diff.num.w$ame95[3] ),
                          cbind(gender = "Woman",outcome = "Difference in Rank Position"  , iv = "Number of \nLand Locations", ame = diff.num.w$ame90[1], ci = 90, lower = diff.num.w$ame90[2], upper = diff.num.w$ame90[3] ),
                          cbind(gender = "Woman",outcome = "Vote Margin for Final Rank"  , iv = "Share of \nKreise Covered", ame = mar.cov.w$ame95[1], ci = 95, lower = mar.cov.w$ame95[2], upper = mar.cov.w$ame95[3] ),
                          cbind(gender = "Woman",outcome = "Vote Margin for Final Rank"  , iv = "Share of \nKreise Covered", ame = mar.cov.w$ame90[1], ci = 90, lower = mar.cov.w$ame90[2], upper = mar.cov.w$ame90[3] ),
                          cbind(gender = "Woman",outcome = "Vote Margin for Final Rank"  , iv = "Number of \nLand Locations", ame = mar.num.w$ame95[1], ci = 95, lower = mar.num.w$ame95[2], upper = mar.num.w$ame95[3] ),
                          cbind(gender = "Woman",outcome = "Vote Margin for Final Rank"  , iv = "Number of \nLand Locations", ame = mar.num.w$ame90[1], ci = 90, lower = mar.num.w$ame90[2], upper = mar.num.w$ame90[3] ),
                        cbind(gender = "Man", outcome = "No. of Ranks Competed"  , iv = "Share of \nKreise Covered", ame = rank.cov.m$ame95[1], ci = 95, lower = rank.cov.m$ame95[2], upper = rank.cov.m$ame95[3] ),
                          cbind(gender = "Man", outcome = "No. of Ranks Competed"  , iv = "Share of \nKreise Covered", ame = rank.cov.m$ame90[1], ci = 90, lower = rank.cov.m$ame90[2], upper = rank.cov.m$ame90[3] ),
                          cbind(gender = "Man", outcome = "No. of Ranks Competed"  , iv = "Number of \nLand Locations", ame = rank.num.m$ame95[1], ci = 95, lower = rank.num.m$ame95[2], upper = rank.num.m$ame95[3] ),
                          cbind(gender = "Man", outcome = "No. of Ranks Competed"  , iv = "Number of \nLand Locations", ame = rank.num.m$ame90[1], ci = 90, lower = rank.num.m$ame90[2], upper = rank.num.m$ame90[3] ),
                          cbind(gender = "Man", outcome = "Difference in Rank Position"  , iv = "Share of \nKreise Covered", ame = diff.cov.m$ame95[1], ci = 95, lower = diff.cov.m$ame95[2], upper = diff.cov.m$ame95[3] ),
                          cbind(gender = "Man", outcome = "Difference in Rank Position"  , iv = "Share of \nKreise Covered", ame = diff.cov.m$ame90[1], ci = 90, lower = diff.cov.m$ame90[2], upper = diff.cov.m$ame90[3] ),
                          cbind(gender = "Man", outcome = "Difference in Rank Position"  , iv = "Number of \nLand Locations", ame = diff.num.m$ame95[1], ci = 95, lower = diff.num.m$ame95[2], upper = diff.num.m$ame95[3] ),
                          cbind(gender = "Man", outcome = "Difference in Rank Position"  , iv = "Number of \nLand Locations", ame = diff.num.m$ame90[1], ci = 90, lower = diff.num.m$ame90[2], upper = diff.num.m$ame90[3] ),
                          cbind(gender = "Man", outcome = "Vote Margin for Final Rank"  , iv = "Share of \nKreise Covered", ame = mar.cov.m$ame95[1], ci = 95, lower = mar.cov.m$ame95[2], upper = mar.cov.m$ame95[3] ),
                          cbind(gender = "Man", outcome = "Vote Margin for Final Rank"  , iv = "Share of \nKreise Covered", ame = mar.cov.m$ame90[1], ci = 90, lower = mar.cov.m$ame90[2], upper = mar.cov.m$ame90[3] ),
                          cbind(gender = "Man", outcome = "Vote Margin for Final Rank"  , iv = "Number of \nLand Locations", ame = mar.num.m$ame95[1], ci = 95, lower = mar.num.m$ame95[2], upper = mar.num.m$ame95[3] ),
                          cbind(gender = "Man", outcome = "Vote Margin for Final Rank"  , iv = "Number of \nLand Locations", ame = mar.num.m$ame90[1], ci = 90, lower = mar.num.m$ame90[2], upper = mar.num.m$ame90[3] )))
                        

row.names(pData) <- 1:nrow(pData)
pData$ame <- as.numeric(as.character(pData$ame))
pData$lower <- as.numeric(as.character(pData$lower))
pData$upper <- as.numeric(as.character(pData$upper))
pData$ci <- as.numeric(as.character(pData$ci))
pData$outcome <- as.character(pData$outcome)
pData$outcome_f <- factor(pData$outcome, levels = c("No. of Ranks Competed","Vote Margin for Final Rank","Difference in Rank Position"))  


library(ggplot2)
#################################################
# Figure 5, Manuscript
#################################################

p = ggplot(pData[pData$ci == 95  ,])
p = p + geom_pointrange(data = pData[pData$ci == 90   ,], 
                        aes(x = iv, y = ame, ymin = lower, ymax = upper),
                        position=position_nudge(x = -0.05), shape = 16)
p = p + geom_pointrange(data = pData[pData$ci == 95   ,], 
                        aes(x = iv, y = ame, ymin = lower, ymax = upper),   lty = 2,
                        position=position_nudge(x = -0.05), shape = 16)
p = p + geom_hline(aes(yintercept = 0), linetype = 2) + xlab("") + ylab("AME")
p = p + facet_grid(gender ~ outcome_f, scales = 'free')
p = p + coord_flip() + theme(axis.text = element_text(size = rel(1), face = "bold", colour = "black"),
                             axis.ticks = element_line(colour = "black"),
                             legend.key = element_rect(fill = "white"),
                             plot.margin = unit(c(0, 0, 0, 0), "cm"),
                             panel.background = element_rect(fill = "white", colour = NA),
                             panel.border = element_rect(fill = NA, colour = "grey50"),
                             strip.background = element_rect(fill = "white", colour = "grey50"),
                             strip.text.x = element_text(size = 12, colour = "black", face = "bold"),
                             strip.text.y = element_text(size = 10, colour = "black", face = "bold"))  
cat("create figure 5\n")
p


#################################################
# explaining local responsiveness
#################################################

 
summary(loc.1 <- lm(pq_both_coverage_land ~ as.factor(land)+ government + rank_before_prob +  age +  office + election + elected + ln_pq_total + margin_final_before + sex
 , data = dat))
loc.1.se.rob <- coeftest(loc.1, vcov = vcovHC(loc.1, "HC"))[,2];coeftest(loc.1, vcov = vcovHC(loc.1, "HC"))
loc.1.se.p <- coeftest(loc.1, vcov = vcovHC(loc.1, "HC"))[,4]
loc.1.se.cl <- coef_test(loc.1, vcov = "CR2", 
                        cluster = loc.1[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,3]
loc.1.se.clp <- coef_test(loc.1, vcov = "CR2", 
                          cluster = loc.1[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,6]

summary(loc.2 <- lm(ln_pq_land ~ as.factor(land)+ government + rank_before_prob + age +  office +  election + elected + ln_pq_total + margin_final_before + sex 
 , data = dat))
loc.2.se.rob <- coeftest(loc.2, vcov = vcovHC(loc.2, "HC"))[,2];coeftest(loc.2, vcov = vcovHC(loc.2, "HC"))
loc.2.se.p <- coeftest(loc.2, vcov = vcovHC(loc.2, "HC"))[,4]
loc.2.se.cl <- coef_test(loc.2, vcov = "CR2", 
                         cluster = loc.2[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,3]
loc.2.se.clp <- coef_test(loc.2, vcov = "CR2", 
                          cluster = loc.2[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,6]

summary(loc.3 <- lm(pq_both_coverage_land ~ as.factor(land)+ government + rank_before_prob +  age +  office + election + elected + ln_pq_total + enec_final_before + sex   
                    , data = dat))
loc.3.se.rob <- coeftest(loc.3, vcov = vcovHC(loc.3, "HC"))[,2];coeftest(loc.3, vcov = vcovHC(loc.3, "HC"))
loc.3.se.p <- coeftest(loc.3, vcov = vcovHC(loc.3, "HC"))[,4]
loc.3.se.cl <- coef_test(loc.3, vcov = "CR2", 
                         cluster = loc.3[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,3]
loc.3.se.clp <- coef_test(loc.3, vcov = "CR2", 
                        cluster = loc.3[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,6]

summary(loc.4 <- lm(ln_pq_land ~ as.factor(land)+ government + rank_before_prob + age +  office +  election + elected + ln_pq_total + enec_final_before + sex 
                    , data = dat))
loc.4.se.rob <- coeftest(loc.4, vcov = vcovHC(loc.4, "HC"))[,2];coeftest(loc.4, vcov = vcovHC(loc.4, "HC"))
loc.4.se.p <- coeftest(loc.4, vcov = vcovHC(loc.4, "HC"))[,4]
loc.4.se.cl <- coef_test(loc.4, vcov = "CR2", 
                         cluster = loc.4[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,3]
loc.4.se.clp <- coef_test(loc.4, vcov = "CR2", 
                        cluster = loc.4[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,6]

summary(loc.5 <- lm(pq_both_coverage_land ~ as.factor(land)+ government + rank_before_prob +  age +  office + election + elected + ln_pq_total + num_cand_before1st + sex   
                    , data = dat))
loc.5.se.rob <- coeftest(loc.5, vcov = vcovHC(loc.5, "HC"))[,2];coeftest(loc.5, vcov = vcovHC(loc.5, "HC"))
loc.5.se.p <- coeftest(loc.5, vcov = vcovHC(loc.5, "HC"))[,4]
loc.5.se.cl <- coef_test(loc.5, vcov = "CR2", 
                         cluster = loc.5[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,3]
loc.5.se.clp <- coef_test(loc.5, vcov = "CR2", 
                        cluster = loc.5[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,6]

summary(loc.6 <- lm(ln_pq_land ~ as.factor(land)+ government + rank_before_prob + age +  office +  election + elected + ln_pq_total + num_cand_before1st + sex 
                    , data = dat))
loc.6.se.rob <- coeftest(loc.6, vcov = vcovHC(loc.6, "HC"))[,2];coeftest(loc.6, vcov = vcovHC(loc.6, "HC"))
loc.6.se.p <- coeftest(loc.6, vcov = vcovHC(loc.6, "HC"))[,4]
loc.6.se.cl <- coef_test(loc.6, vcov = "CR2", 
                         cluster = loc.6[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,3]
loc.6.se.clp <- coef_test(loc.6, vcov = "CR2", 
                     cluster = loc.6[["model"]]$`as.factor(land)`, test = "Satterthwaite")[,6]


#################################################
# Table 12 in Appendix
#################################################
library(texreg)
cat("create table 12\n")
texreg(list(loc.1, loc.1, loc.2, loc.2,loc.3, loc.3, loc.4, loc.4,loc.5, loc.5, loc.6, loc.6),
       override.se = list(loc.1.se.rob, loc.1.se.cl, loc.2.se.rob, loc.2.se.cl,
                    loc.3.se.rob, loc.3.se.cl, loc.4.se.rob, loc.4.se.cl,
                    loc.5.se.rob, loc.5.se.cl, loc.6.se.rob, loc.6.se.cl), 
       override.pvalues = list(loc.1.se.p, loc.1.se.clp, loc.2.se.p, loc.2.se.clp,
                               loc.3.se.p, loc.3.se.clp, loc.4.se.p, loc.4.se.clp,
                               loc.5.se.p, loc.5.se.clp, loc.6.se.p, loc.6.se.clp), 
       omit.coef = "land", 
       reorder.coef = c(2,3,4,5,6,7,8,9,11,12,10,1),
          custom.coef.names = c( "Constant","Government", "Pr (Election) Rank t-1",
                                "Age", "Office Holder","Time Trend (Election Year)",
                                "Times Elected", "Total No. of PQs, log","Margin on Final Rank t-1", "Man",
                                "Eff. No. of Cand. for Final Rank t-1",
                                "No. of Candidates for Final Rank t-1"),
          stars = c(0.1, 0.05, 0.01), digits = 3)



#################################################
# 1) number of ranks competed
#################################################
##### -- Women -- #####
### --- PQs - Number --- ###
summary(w2.d <- glm(ranks_competed ~ -1+as.factor(incumbent)+  rank_before_prob + 
                      elected+ theta_mean_land_deviation + ln_pq_total  + ln_pq_other,  data = dat[dat$sex == "female",], family = poisson("log")))
w2.d.se.rob <- coeftest(w2.d, vcov = vcovHC(w2.d, "HC"))[,2];coeftest(w2.d, vcov = vcovHC(w2.d, "HC"))
w2.d.se.cl <- coef_test(w2.d, vcov = "CR2",  cluster = dat$incumbent[dat$sex == "female"], test = "Satterthwaite")[,3]

##### -- Men -- #####
### --- PQs - Number --- ###
summary(m2.d <- glm(ranks_competed ~ -1+as.factor(incumbent) + rank_before_prob + 
                      office + elected+ theta_mean_land_deviation +ln_pq_total +  ln_pq_other,  data = dat[dat$sex == "male",], family = poisson("log")))
m2.d.se.rob <- coeftest(m2.d, vcov = vcovHC(m2.d, "HC"))[,2];coeftest(m2.d, vcov = vcovHC(m2.d, "HC"))
m2.d.se.cl <- coef_test(m2.d, vcov = "CR2",  cluster = dat$incumbent[dat$sex == "male"], test = "Satterthwaite")[,3]

#################################################
# 2) difference in rank
#################################################

##### -- Women -- #####
### --- PQs - Number --- ###
summary(w4.d <- lm(rank_diff ~ -1+as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation + ln_pq_total +  ln_pq_other,  data = dat[dat$sex == "female",]))
w4.d.se.rob <- coeftest(w4.d, vcov = vcovHC(w4.d, "HC"))[,2];coeftest(w4.d, vcov = vcovHC(w4.d, "HC"))
w4.d.se.cl <- coef_test(w4.d, vcov = "CR2", 
                        cluster = dat$incumbent[dat$sex == "female"], test = "Satterthwaite")[,3]

##### -- Men -- #####
### --- PQs - Number --- ###
summary(m4.d <- lm(rank_diff ~ -1+as.factor(incumbent) + rank_before_prob + 
                     office + elected+ theta_mean_land_deviation+ ln_pq_total +  ln_pq_other,  data = dat[dat$sex == "male",]))#, family = binomial("logit")))
m4.d.se.rob <- coeftest(m4.d, vcov = vcovHC(m4.d, "HC"))[,2];coeftest(m4.d, vcov = vcovHC(m4.d, "HC"))
m4.d.se.cl <- coef_test(m4.d, vcov = "CR2",    cluster = dat$incumbent[dat$sex == "male"], test = "Satterthwaite")[,3]


#################################################
# 3) margin for final rank
#################################################

##### -- Women -- #####
### --- PQs - Number --- ###
summary(w1.d <- lm(margin_final ~ -1+as.factor(incumbent) + rank_before_prob +
                     office + elected+ theta_mean_land_deviation + ln_pq_total +  ln_pq_other,  data = dat[dat$sex == "female",]))
w1.d.se.rob <- coeftest(w1.d, vcov = vcovHC(w1.d, "HC"))[,2];coeftest(w1.d, vcov = vcovHC(w1.d, "HC"))
w1.d.se.cl <- coef_test(w1.d, vcov = "CR2",    cluster = dat$incumbent[dat$sex == "female"], test = "Satterthwaite")[,3]


##### -- Men -- #####
### --- PQs - Number --- ###
summary(m1.d <- lm(margin_final ~ -1+as.factor(incumbent) + rank_before_prob +
                     office + elected+ theta_mean_land_deviation + ln_pq_total +  ln_pq_other,  data = dat[dat$sex == "male",]))
m1.d.se.cl <- coef_test(m1.d, vcov = "CR2",  cluster = dat$incumbent[dat$sex == "male"], test = "Satterthwaite")[,3]
m1.d.se.rob <- coeftest(m1.d, vcov = vcovHC(m1.d, "HC"))[,2];coeftest(m1.d, vcov = vcovHC(m1.d, "HC"))



#################################################
# Table 14 in the Appendix
#################################################
# split by gender
cat("create table 14\n")
stargazer(list(  w2.d, m2.d, w4.d, m4.d,w1.d, m1.d
                 ),
          se = list( w2.d.se.rob, m2.d.se.rob, w4.d.se.rob, m4.d.se.rob,w1.d.se.rob, m1.d.se.rob
                     ), omit = "incumbent",
          covariate.labels = c("Pr (Election) Rank t-1",
                               "Office Holder",
                               "Times Elected",
                               "Ideological Deviation from Landesverband",
                               "Total Number of PQs, log",
                               "Number of Locations outside of PR District Mentioned in PQs, log"))

#################################################
# Table 15 in the Appendix
#################################################
# split by gender
cat("create table 15\n")
stargazer(list( 
                 w2.d, m2.d, w4.d, m4.d,w1.d, m1.d),
          se = list( 
                     w2.d.se.cl, m2.d.se.cl, w4.d.se.cl, m4.d.se.cl,w1.d.se.cl, m1.d.se.cl), 
          omit = "incumbent",
          covariate.labels = c("Pr (Election) Rank t-1",
                               "Office Holder",
                               "Times Elected",
                               "Ideological Deviation from Landesverband",
                               "Total Number of PQs, log",
                               "Number of Locations outside of PR District Mentioned in PQs, log"))





#################################################
# Distribution of PQs
#################################################
#################################################
# Figure 13 in the Appendix
#################################################
cat("create figure 13\n")
par(mfrow = c(2,1),
    oma = c(1.5,1,1,0) + 0.1,
    mar = c(1.5,1,2,1) + 0.1)
hist(dat$pq_total, breaks = 50,  xlab = "", main = "Total No. of PQs", xaxt = "n", yaxt = "n", col = "white")
axis(1, seq(0,750,100), seq(0,750,100), col.ticks = "black", col = "white")
axis(2, seq(0,180,20), seq(0,180,20), col.ticks = "black", col = "white")

hist(dat$pq_both_num_land, breaks = 50,  xlab = "", main = "No. of Land Locations Mentioned in PQs", xaxt = "n", yaxt = "n", col = "white")
axis(1, seq(0,250,50), seq(0,250,50), col.ticks = "black", col = "white")
axis(2, seq(0,200,20), seq(0,200,20), col.ticks = "black", col = "white")

cat("end logging")
sink()
